home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / ADI.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  4KB  |  128 lines
  1. PROCEDURE adi(a,b,c,d,e,f,g: gljmax; VAR u: gljmax;
  2.       jmax,k: integer; alpha,beta,eps: double);
  3. (* Programs using routine ADI must define the type
  4. TYPE
  5.    gljmax = ARRAY [1..jmax,1..jmax] OF double;
  6. in the main routine. *)
  7. LABEL 99;
  8. CONST
  9.    jj=50;
  10.    kk=6;
  11.    nrr=32;      (* nrr=2 to the power (kk-1) *)
  12.    maxits=100;
  13.    zero=0.0;
  14.    two=2.0;
  15.    half=0.5;
  16. TYPE
  17.    gljjarray = ARRAY [1..jj] OF double;
  18. VAR
  19.    i,nr,nits,next,n,l,kits,k1,j,twopwr: integer;
  20.    rfact,resid,disc,anormg,anorm,ab: double;
  21.    aa,bb,cc,rr,uu: gljjarray;
  22.    psi: ARRAY [1..jj,1..jj] OF double;
  23.    alph,bet: ARRAY [1..kk] OF double;
  24.    r: ARRAY [1..nrr] OF double;
  25.    s: ARRAY [1..nrr,1..kk] OF double;
  26. PROCEDURE tridag(a,b,c,r: gljjarray; VAR u: gljjarray; n: integer);
  27.    (* This is a double precision version of TRIDAG for use with ADI,
  28.       which defines the constant jj and the type gljjarray.   *)
  29.    VAR
  30.       j: integer;
  31.       bet: double;
  32.       gam: gljjarray;
  33.    BEGIN
  34.       IF (b[1] = 0.0) THEN BEGIN
  35.          writeln('pause 1 in TRIDAG'); readln END;
  36.       bet := b[1];
  37.       u[1] := r[1]/bet;
  38.       FOR j := 2 TO n DO BEGIN
  39.          gam[j] := c[j-1]/bet;
  40.          bet := b[j]-a[j]*gam[j];
  41.          IF (bet = 0.0) THEN BEGIN
  42.             writeln('pause 2 in TRIDAG'); readln END;
  43.          u[j] := (r[j]-a[j]*u[j-1])/bet
  44.       END;
  45.       FOR j := n-1 DOWNTO 1 DO u[j] := u[j]-gam[j+1]*u[j+1]
  46.    END;
  47. BEGIN
  48.    IF (jmax > jj) THEN BEGIN
  49.       writeln('Pause in routine ADI - increase jj'); readln
  50.    END;
  51.    IF (k > (kk-1)) THEN BEGIN
  52.       writeln('Pause in routine ADI - increase kk'); readln
  53.    END;
  54.    k1 := k+1;
  55.    nr := 1;
  56.    FOR i := 1 TO k DO nr := 2*nr;
  57.    alph[1] := alpha;
  58.    bet[1] := beta;
  59.    FOR j := 1 TO k DO BEGIN
  60.       alph[j+1] := sqrt(alph[j]*bet[j]);
  61.       bet[j+1] := half*(alph[j]+bet[j])
  62.    END;
  63.    s[1,1] := sqrt(alph[k1]*bet[k1]);
  64.    FOR j := 1 TO k DO BEGIN
  65.       ab := alph[k1-j]*bet[k1-j];
  66.       twopwr := 1;
  67.       FOR i := 1 TO (j-1) DO twopwr := 2*twopwr;
  68.       FOR n := 1 TO twopwr DO BEGIN
  69.          disc := sqrt(sqr(s[n,j])-ab);
  70.          s[2*n,j+1] := s[n,j]+disc;
  71.          s[2*n-1,j+1] := ab/s[2*n,j+1]
  72.       END
  73.    END;
  74.    FOR n := 1 TO nr DO r[n] := s[n,k1];
  75.    anormg := zero;
  76.    FOR j := 2 TO (jmax-1) DO BEGIN
  77.       FOR l := 2 TO (jmax-1) DO BEGIN
  78.          anormg := anormg+abs(g[j,l]);
  79.          psi[j,l] := -d[j,l]*u[j,l-1]
  80.             +(r[1]-e[j,l])*u[j,l]-f[j,l]*u[j,l+1]
  81.       END
  82.    END;
  83.    nits := maxits DIV nr;
  84.    FOR kits := 1 TO nits DO BEGIN
  85.       FOR n := 1 TO nr DO BEGIN
  86.          IF (n = nr) THEN next := 1 ELSE next := n+1;
  87.          rfact := r[n]+r[next];
  88.          FOR l := 2 TO (jmax-1) DO BEGIN
  89.             FOR j := 2 TO jmax-1 DO BEGIN
  90.                aa[j-1] := a[j,l];
  91.                bb[j-1] := b[j,l]+r[n];
  92.                cc[j-1] := c[j,l];
  93.                rr[j-1] := psi[j,l]-g[j,l]
  94.             END;
  95.             tridag(aa,bb,cc,rr,uu,jmax-2);
  96.             FOR j := 2 TO (jmax-1) DO BEGIN
  97.                psi[j,l] := -psi[j,l]
  98.                   +two*r[n]*uu[j-1]
  99.             END
  100.          END;
  101.          FOR j := 2 TO (jmax-1) DO BEGIN
  102.             FOR l := 2 TO (jmax-1) DO BEGIN
  103.                aa[l-1] := d[j,l];
  104.                bb[l-1] := e[j,l]+r[n];
  105.                cc[l-1] := f[j,l];
  106.                rr[l-1] := psi[j,l]
  107.             END;
  108.             tridag(aa,bb,cc,rr,uu,jmax-2);
  109.             FOR l := 2 TO (jmax-1) DO BEGIN
  110.                u[j,l] := uu[l-1];
  111.                psi[j,l] := -psi[j,l]+rfact*uu[l-1]
  112.             END
  113.          END
  114.       END;
  115.       anorm := zero;
  116.       FOR j := 2 TO (jmax-1) DO BEGIN
  117.          FOR l := 2 TO (jmax-1) DO BEGIN
  118.             resid := a[j,l]*u[j-1,l]+(b[j,l]+e[j,l])*u[j,l]
  119.                +c[j,l]*u[j+1,l]+d[j,l]*u[j,l-1]
  120.                +f[j,l]*u[j,l+1]+g[j,l];
  121.             anorm := anorm+abs(resid)
  122.          END
  123.       END;
  124.       IF (anorm < (eps*anormg)) THEN GOTO 99
  125.    END;
  126.    writeln('Pause in routine ADI - too many iterations');
  127. 99:   END;
  128.